Sanjeeva Reddy Dodlapati

LinkedIn GitHub Google Scholar Twitter Blog Email 
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import dash_bio
import scipy
import math
data = pd.read_parquet('max_gve.parquet.gzip')
data
chrom pos name abs_max_effect max_column ave_effect
0 1 174946171 rs147831690 0.027101 nrnlsc_H3K79me1 0.000700
1 1 174946180 rs1389168246 0.024842 nrnlsc_H3K79me1 0.000408
2 1 174946196 rs1374951931 0.072019 nrnlsc_H3K79me1 0.001219
3 1 174946202 rs948828419 0.033973 nrnlsc_H3K79me1 0.001874
4 1 174946211 rs896018724 0.064587 nrnlsc_H3K79me1 0.001534
... ... ... ... ... ... ...
547681 2 156681080 rs747765225 0.009990 SKNSH_ARNTL 0.000462
547682 2 156681082 rs994964368 0.009063 cdtn_H3K4me3 0.000822
547683 2 156681084 rs1028224262 0.008285 SKNSH_ARNTL 0.000428
547684 2 156681085 rs541069773 0.016144 SKNSH_ARNTL 0.000941
547685 2 156681113 rs1432010733 0.016751 MOC_NP_ATAC 0.000612

16476152 rows × 6 columns

dfs = []
for i in range(10):
    print(i)
    df = data.iloc[:1000, :]
    dfs.append(df)
print(df.head())
0
1
2
3
4
5
6
7
8
9
   chrom        pos          name  abs_max_effect       max_column  ave_effect
0      1  174946171   rs147831690        0.027101  nrnlsc_H3K79me1    0.000700
1      1  174946180  rs1389168246        0.024842  nrnlsc_H3K79me1    0.000408
2      1  174946196  rs1374951931        0.072019  nrnlsc_H3K79me1    0.001219
3      1  174946202   rs948828419        0.033973  nrnlsc_H3K79me1    0.001874
4      1  174946211   rs896018724        0.064587  nrnlsc_H3K79me1    0.001534
dfs = pd.concat(dfs)
dfs
chrom pos name abs_max_effect max_column ave_effect
0 1 174946171 rs147831690 0.027101 nrnlsc_H3K79me1 0.000700
1 1 174946180 rs1389168246 0.024842 nrnlsc_H3K79me1 0.000408
2 1 174946196 rs1374951931 0.072019 nrnlsc_H3K79me1 0.001219
3 1 174946202 rs948828419 0.033973 nrnlsc_H3K79me1 0.001874
4 1 174946211 rs896018724 0.064587 nrnlsc_H3K79me1 0.001534
... ... ... ... ... ... ...
995 1 174949547 rs1573971060 0.042717 MOC_NP_ATAC 0.001534
996 1 174949551 rs1158195944 0.005916 SKNSH_H3F3A 0.000336
997 1 174949556 rs191506112 0.006997 ACC_NP_ATAC 0.000302
998 1 174949557 rs557939733 0.005723 ACC_NP_ATAC 0.000226
999 1 174949559 rs551299464 0.007402 MiddleFrGy_DNase 0.000275

10000 rows × 6 columns

data = pd.read_csv('max_gve.csv')
data['chrom'][data.chrom=='X'] = 23
data['chrom'][data.chrom=='Y'] = 24
data['chrom'][data.chrom=='16.0'] = 16
data['chrom'] = data.chrom.astype(int)
data
chrom pos name max_effect max_column ave_effect
0 1 238807165 rs1207962803 0.015620 brnmec_DNase 0.000920
1 1 238807165 rs1242649293 0.081236 nrnlsc_H3K9ac -0.002660
2 1 238807167 rs997302978 0.002643 ACC_NP_ATAC -0.001335
3 1 238807168 rs1270214146 0.045428 nrnlsc_H3K79me1 -0.000773
4 1 238807168 rs1270214146 0.083614 nrnlsc_H3K79me1 0.001402
... ... ... ... ... ... ...
16652518 2 156681080 rs747765225 0.006044 PMC_NN_ATAC -0.000144
16652519 2 156681082 rs994964368 0.000293 MeT_NN_ATAC -0.000791
16652520 2 156681084 rs1028224262 0.004548 VPC_NN_ATAC -0.000109
16652521 2 156681085 rs541069773 0.011259 VPC_NN_ATAC -0.000122
16652522 2 156681113 rs1432010733 0.010121 nrnlsc_H3K4me1 0.000242

16652523 rows × 6 columns

data.chrom.unique()
array([ 1, 22,  3, 12, 16, 17, 13, 10, 11,  2, 19, 14, 15, 21])
neg = -(data.max_effect.values)
exp = np.exp(neg)
data['exp'] = exp

zscores = scipy.stats.zscore(data.ave_effect.values)
data['zscore'] = zscores
pvalues = scipy.stats.norm.sf(zscores)
data['p'] = pvalues

data
chrom pos name max_effect max_column ave_effect exp zscore p
0 1 238807165 rs1207962803 0.015620 brnmec_DNase 0.000920 0.984502 0.310131 0.378231
1 1 238807165 rs1242649293 0.081236 nrnlsc_H3K9ac -0.002660 0.921976 -0.931275 0.824144
2 1 238807167 rs997302978 0.002643 ACC_NP_ATAC -0.001335 0.997360 -0.471795 0.681464
3 1 238807168 rs1270214146 0.045428 nrnlsc_H3K79me1 -0.000773 0.955588 -0.276915 0.609077
4 1 238807168 rs1270214146 0.083614 nrnlsc_H3K79me1 0.001402 0.919786 0.477353 0.316555
... ... ... ... ... ... ... ... ... ...
16652518 2 156681080 rs747765225 0.006044 PMC_NN_ATAC -0.000144 0.993974 -0.058674 0.523394
16652519 2 156681082 rs994964368 0.000293 MeT_NN_ATAC -0.000791 0.999707 -0.283229 0.611499
16652520 2 156681084 rs1028224262 0.004548 VPC_NN_ATAC -0.000109 0.995462 -0.046438 0.518520
16652521 2 156681085 rs541069773 0.011259 VPC_NN_ATAC -0.000122 0.988804 -0.050923 0.520306
16652522 2 156681113 rs1432010733 0.010121 nrnlsc_H3K4me1 0.000242 0.989930 0.075191 0.470031

16652523 rows × 9 columns

data.to_csv('abs_max_ave_effect_pvalue.csv', index=False)
data = pd.read_csv('tbln_H3K27ac_kmeans0.csv')
data['chrom'][data.chrom=='X'] = 23
data['chrom'][data.chrom=='Y'] = 24
data['chrom'][data.chrom=='16.0'] = 16
data['chrom'] = data.chrom.astype(int)
data
chrom pos name tbln_H3K27ac diff abs_diff
0 10 84139514.0 rs932821133 0.105500 0.105500 0.105500
1 10 84139543.0 rs1229596357 -0.183456 -0.183456 0.183456
2 10 84139573.0 rs1481839079 0.143182 0.143182 0.143182
3 10 84140064.0 rs1452183496 0.136644 0.136644 0.136644
4 10 84140294.0 rs1429667538 0.100358 0.100358 0.100358
... ... ... ... ... ... ...
726226 3 194089744.0 rs1351149363 -0.218239 -0.218239 0.218239
726227 3 194089846.0 rs890061713 -0.116548 -0.116548 0.116548
726228 3 194089911.0 rs1246306491 -0.164336 -0.164336 0.164336
726229 3 194090129.0 rs1181593567 -0.111510 -0.111510 0.111510
726230 3 194095916.0 rs1009672135 -0.116104 -0.116104 0.116104

726231 rows × 6 columns

neg = -(data.tbln_H3K27ac.values)
exp = np.exp(neg)
data['exp'] = exp

zscores = scipy.stats.zscore(data.tbln_H3K27ac.values)
data['zscore'] = zscores
pvalues = scipy.stats.norm.sf(zscores)
data['p'] = pvalues

data
chrom pos name tbln_H3K27ac diff abs_diff exp
0 10 84139514.0 rs932821133 0.105500 0.105500 0.105500 0.899874
1 10 84139543.0 rs1229596357 -0.183456 -0.183456 0.183456 1.201363
2 10 84139573.0 rs1481839079 0.143182 0.143182 0.143182 0.866597
3 10 84140064.0 rs1452183496 0.136644 0.136644 0.136644 0.872280
4 10 84140294.0 rs1429667538 0.100358 0.100358 0.100358 0.904514
... ... ... ... ... ... ... ...
726226 3 194089744.0 rs1351149363 -0.218239 -0.218239 0.218239 1.243884
726227 3 194089846.0 rs890061713 -0.116548 -0.116548 0.116548 1.123611
726228 3 194089911.0 rs1246306491 -0.164336 -0.164336 0.164336 1.178611
726229 3 194090129.0 rs1181593567 -0.111510 -0.111510 0.111510 1.117965
726230 3 194095916.0 rs1009672135 -0.116104 -0.116104 0.116104 1.123113

726231 rows × 7 columns

plt.hist(data.tbln_H3K27ac)
(array([2.40000e+01, 5.00000e+02, 1.43580e+04, 3.43638e+05, 0.00000e+00,
        3.38365e+05, 2.81640e+04, 1.12700e+03, 4.00000e+01, 1.50000e+01]),
 array([-0.74844134, -0.58586169, -0.42328205, -0.2607024 , -0.09812275,
         0.0644569 ,  0.22703654,  0.38961619,  0.55219584,  0.71477548,
         0.87735513]),
 <BarContainer object of 10 artists>)

# scipy.stats.t.sf(data.ITC_NN_ATAC.values, df=2)
# scipy.stats.t.sf(data.ITC_NN_ATAC.values, df=3)
chrom pos name tbln_H3K27ac diff abs_diff exp zscore p
0 10 84139514.0 rs932821133 0.105500 0.105500 0.105500 0.899874 0.658540 0.255096
1 10 84139543.0 rs1229596357 -0.183456 -0.183456 0.183456 1.201363 -1.181742 0.881346
2 10 84139573.0 rs1481839079 0.143182 0.143182 0.143182 0.866597 0.898519 0.184454
3 10 84140064.0 rs1452183496 0.136644 0.136644 0.136644 0.872280 0.856886 0.195754
4 10 84140294.0 rs1429667538 0.100358 0.100358 0.100358 0.904514 0.625787 0.265727
... ... ... ... ... ... ... ... ... ...
726226 3 194089744.0 rs1351149363 -0.218239 -0.218239 0.218239 1.243884 -1.403260 0.919730
726227 3 194089846.0 rs890061713 -0.116548 -0.116548 0.116548 1.123611 -0.755621 0.775062
726228 3 194089911.0 rs1246306491 -0.164336 -0.164336 0.164336 1.178611 -1.059971 0.855421
726229 3 194090129.0 rs1181593567 -0.111510 -0.111510 0.111510 1.117965 -0.723538 0.765325
726230 3 194095916.0 rs1009672135 -0.116104 -0.116104 0.116104 1.123113 -0.752797 0.774214

726231 rows × 9 columns

data.chrom.unique()
array([10,  1, 16, 23, 15, 21,  2,  6,  7,  5, 12, 19,  8,  3,  9,  4, 11,
       24, 13, 14, 17, 22, 20, 18])
data.to_csv('data_with_pvalue.csv', index=False)
dash_bio.ManhattanPlot(dataframe=df, chrm='CHR', bp='BP', snp='SNP', p='P')
file_lists = [files[:500], files[500:1000], files[1000:1500], files[1500:2000], files[2000:2500], files[2500:3000],
              files[3000:3500], files[3500:4000], files[4000:4500], files[4500:5000], files[5000:]]
df = pd.read_csv('https://raw.githubusercontent.com/plotly/dash-bio-docs-files/master/manhattan_data.csv')
df
CHR BP P SNP ZSCORE EFFECTSIZE GENE DISTANCE
0 1 937641 0.335344 rs9697358 0.9634 -0.0946 ISG15 1068
1 1 1136887 0.245857 rs34945898 1.1605 -0.0947 TNFRSF4 0
2 1 2116240 0.823286 rs12034613 0.2233 -0.0741 FP7162 0
3 1 2310562 0.493204 rs4648633 0.6852 0.0146 MORN1 0
4 1 2681715 0.605392 rs4430271 0.5167 0.1234 MMEL1 127427
... ... ... ... ... ... ... ... ...
14407 23 153207974 0.574920 rs766418 0.5608 -0.0190 NaN -1
14408 23 153280339 0.978400 rs11593 0.0271 -0.1355 NaN -1
14409 23 153546061 0.056012 rs5987005 1.9109 0.0026 NaN -1
14410 23 153903746 0.132189 rs28370194 1.5055 -0.1425 NaN -1
14411 23 154208789 0.497909 rs524609 0.6778 0.2542 NaN -1

14412 rows × 8 columns

# dash_bio.ManhattanPlot(dataframe=data, chrm='CHR', bp='pos')
files = [i for i in range(5561)]
len(files)
5561
i = 0
step = 500
file_lists = []
while i+step <= len(files):
    print(i)
    flist = files[i:i+step]
    i+=step
    file_lists.append(flist)
    
flist = files[i:]
file_lists.append(flist)
print(len(file_lists))
0
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
12
len(file_lists[-1])
61